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By simulating the first order globule-crystal transition of a flexible homopolymer chain, both 
by collision dynamics and Monte Carlo with non-kinetic moves, we show that the effective and the 
thermodynamic transition temperatures are different and we propose a way of quantifying the kinetic 
hindering. We then also observe that the top eigenvalue in the spectrum of the dynamical (contact 
or adjacency) matrix provides insight into the ensembles of folding and unfolding trajectories, and 
may be a suitable additional reaction coordinate for the folding transition of chain molecules. 
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Dynamical simulations of phase transitions in simple 
chain molecules are important for understanding the un- 
derlying microscopic folding mechanisms of polymers and 
proteins. These simulations are usually computationally 
expensive, because the phases are often separated by a 
large free energy barrier, and the fluctuation leading to 
barrier crossing is a rare event on the molecular simula- 
tion timescale. Free energy barriers are usually studied 
by Monte Carlo (MC) simulation, and dynamical stud- 
ies of barrier crossing are less common. Here, we use 
stochastic hard-sphere molecular dynamics, collision dy- 
namics (CD) for short, to study the kinetics involved 
in the first-order crystallisation transition for a flexible 
homopolymer model, composed of bonded hard spheres 
with square-well nonbonded attractions, which has re- 
cently been intensively studied by MC [l|-|5|. By using 
CD, we avoid unphysical aspects of Monte Carlo (such 
as connectivity-altering moves improving sampling effi- 
ciency) and hence are able to observe the effects of ki- 
netic bottlenecks in determining the rate, for a situation 
that mimics realistic dynamics. 

In what follows, the polymer is represented by 128 
hard spheres with two types of square-well potential for 
bonded and non-bonded pairs of monomers: 
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where is the distance between centers of two 
monomers i and j, a is the diameter of the bead, e is 
the depth of the square well, and x an d Xb are the rela- 
tive widths of the square well and the nearest-neighbour 
bond respectively. In the MC studies 0] Xb = 1, but 
for simple hard sphere dynamics, we choose a slightly 
larger value Xb = 1-04; we have verified that this makes 
very little difference to the equilibrium properties. Taylor 
et al. 0, [H have determined that an all-in-one 'protein- 
like' crystallisation from the expanded state occurs for 



X < 1-06, and a two-step 'polymer-like' mechanism via 
a liquid-like globule, for x ^ 1.06; we study the globule- 
crystal stage of the two-step process for x values in this 
vicinity. All beads have equal mass m which we take 
equal to unity. Throughout, we work in reduced units: 
fj = 1, e = 1 and fee = 1 (Boltzmann's constant) so 
T = k^T/e. Corresponding real values, for monomer 
beads corresponding to amino acids of the kind found in 
proteins, would be m « 2 x 10~ 25 kg, a « 6 x 10~ 10 m, 
e w 7 x 10~ 22 J, and a time unit rj 1CP 11 s. 



We have simulated the above chain model with both 
the WL and CD methods. The WL Monte Carlo move 
set @, 0] consists of crankshaft, pivot, end-bridging, and 
regrowth moves; the latter two being connectivity alter- 
ing. The regrowth move consists of regrowing up to 3 
beads at either end of the chain, using a configuration- 
bias algorithm, and includes the possibility of reversing 
the chain. This MC move set was used with the WL al- 
gorithm as in Q to iteratively approximate the density- 
of-states function g(E), giving a well-sampled set of con- 
figurations across the whole energy range. In CD, free 
flight of the spheres occurs between elastic collisions in 
standard fashion 0, [j| . Collisions occur at each discon- 
tinuity in cqns ([2]). Additionally, thermal jolts res- 
elect the velocities of individual atoms from the Boltz- 
mann distribution and introduce a stochasticity into the 
dynamics. The time separation of the jolts has a Pois- 
son distribution with mean time r giving the strength of 
the coupling (typically r = 0.1 in reduced units). The 
WL simulations were used to determine the thermody- 
namic freezing temperatures, Tf, at which doubly-peaked 
canonical ensemble energy distributions Pc{E) were ob- 
served, for a range of x values. A typical example is 
shown in Figure [TJ Similar probability distributions were 
determined from the CD simulations. These confirmed a 
unique globule phase; we denote this as B and identify 
a single energy distribution Pg(E). However, CD sim- 
ulations initialized in different realisations of the crystal 
phase did not explore as wide an energy range as the 
corresponding WL simulations. From this we infer that 
the crystal phase (at Tf) consists of a large number of 
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FIG. 1: Free energy F c (-B) = E - T\n{g(E)) (black) and 
canonical probability function Pc(E) oc g(E)exp(—E/kBT) 
(green) obtained from WL simulation. Vertical lines schema- 
tize the interfaces Xf in folding (red) and the interfaces Xf 
in unfolding (blue) directions. Displayed also are energy dis- 
tributions of states in the crystalline phase (blue) and in the 
globule phase (red); these graphs must be scaled according 
to the folding and unfolding rates [l(| to be comparable with 
Pc(E). The vertical axis of the probability function is not 
shown; these functions are normalized. The data are for the 
chain with X = 1-07, T = 0.498. 



FIG. 2: Chevron plots with folding and unfolding rates com- 
puted by FFS. The intersections give the transition tempera- 
tures Tf u estimated by collision dynamics. The schematic for 
X = 1-06 shows that kinetic hindering of unfolding in CD sim- 
ulations can explain the observation that Tf D > Tf 10 . The 
schematic assumes that the hindering involved in the folding 
process is negligible. 



choice of order parameter. The rate from crystal to glob- 
ule is given by: 



basins with slightly different mean energies, separated 
by kinetic barriers which cannot be overcome (at these 
temperatures) on a simulation timescalc without the use 
of unphysical MC moves. We therefore identify the crys- 
tal state A as being the state, that can be reached from 
other states within the distribution function Pa(E) via 
constant temperature MC (including unphysical moves), 
which has the lowest mean energy. CD alone succeeded 
in sampling these low-lying basins only for x _! 1-05. 

The order parameter A describing the qualitative dif- 
ference between the globule and crystal state was cho- 
sen to be equal to the potential energy E of the chain. 
To accelerate the sampling of fluctuations leading to the 
folding and unfolding transitions we used Forward Flux 
Sampling (FFS) [TT - 13 , which separates the phase space 
by n hypcrplancs orthogonal to the order parameter A 
and measures the probability flux through these planes. 
Here, the hyperplanes (also called interfaces) and the as- 
sociated energies are denoted by the same symbol Aj. 
FFS for our chain was performed in both directions, 
i.e. from globule to crystal, and from crystal to globule. 
This approach allows us to focus on kinetic effects associ- 
ated with each direction separately, and is better suited 
to our dynamics than the backwards/forwards shooting 
approaches of alternatives such as Transition Interface 
Sampling, which are otherwise essentially equivalent [13 ]. 
FFS is also known to be relatively insensitive to the 
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Here x aj is the probability flux through Aq ; state 

A is defined by an interface Xa < Aq 4 . (Ha) is the prob- 
ability of being most recently in state A. The fraction 

0&A I (^a) then simply represents the inverse of the 

average time needed to reach Xq from the first crossing of 
A^. P{Xf +1 | Xf) is the conditional probability of reach- 
ing the interface Xf +1 from Xf , which is given by the 
fraction of partial pathways started from Xf which reach 
Xf +1 before they fall back to A^. The partial pathways 
are eventually connected into the full transition pathways 
starting in A and ending in B. The rate ks^A from glob- 
ule to crystal is given analogously. In this way, we gath- 
ered 8192 folding and unfolding transition pathways by 
CD simulation for each x an d T, choosing temperatures 
in the vicinity of the transition temperatures determined 
by WL. 

The interfaces in FFS are positioned as follows. Ener- 
gies Aq and Aq 3 are chosen such that the ranges A < Aq 
and X > Xq capture 99.9% of the corresponding inte- 
grated densities Pa and Pb ■ The boundaries of A and B 
are then defined as Xa = Xq— 50 and A# = Aq 3 + 50. To 
reconstruct the full pathways wc define an extra plane as 
Aj 4 = Aq 3 . Interface X^_ l is placed close to the isocom- 
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FIG. 3: Phase diagram with freezing temperatures deter 

iMC 



mined from chevron plots {Tf° , FFS and BF) and by 



Maxwell construction (Tf , WL). The errors are smaller than 
the symbol sizes. 



mittor, specifically such that P(X B \ X^-i) = 0.9. The 
remaining interfaces X^, are chosen according to 
the optimization scheme of Ref . [HI . An analogous pro- 
cedure applies to interfaces Af , A^. 

To avoid undersampling of X B , 32 random disordered 
chains are equilibrated, and X B is sampled in parallel 
starting from these different chains. Similarly, CD simu- 
lations started from different points in A^ gave P < f > {E) 1 
i = 1, 256. Aq is then sampled in parallel starting from 
a configuration belonging to the distribution Pa (E) hav- 
ing the lowest mean energy, and from 31 other configura- 
tions which are separated from this state by a large num- 
ber of MC moves, including connectivity-altering moves, 
conducted at constant temperature. 

Transition temperatures Tp D , defined dynamically by 
the equation Ua^b = kB^A, are obtained from the 
FFS simulations via a so-called chevron plot, Fig. [5] 
(based on the assumption of Arrhcnius-likc behaviour). 
These were slightly, but systematically, higher than the 
WL transition temperatures T^ 10 by an amount AT 
t cd _ T ' 
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0.005-0.012 (Fig.©. For the larger values 
of x, barrier crossing could be observed directly, and rate 
constants were calculated by brute force (BF) simulation. 
Transition temperatures obtained in this way were again 
comparable to the temperatures obtained in MC simu- 
lation, but were also systematically higher. This sug- 
gests that the discrepancy is due to a real dynamical (ki- 
netic) effect rather than any deficiency of the FFS algo- 
rithm itself. The most likely explanation is that CD, hin- 
dered by metastable basins acting as kinetic traps, yields 
significantly lower unfolding rates than the MC simula- 
tion, which escapes the traps with the aid of non-kinetic 
connectivity-altering moves. Fig. [2] indicates schemati- 



cally the extent to which the unfolding regression line 
must be shifted to give the observed shift AT; this could 
provide a quantitative measure of the kinetic hindering. 

The remaining part of the paper points out a remark- 
able property the dynamical matrix of our polymer sys- 
tem in the vicinity of the free energy maximum. We focus 
on the chain with \ = 1-07 and T = 0.498, but similar re- 
sults were obtained for other parameters, becoming even 
more distinct with decreasing \- The dynamical matrix 
of our chain is defined as follows: 
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if i ^ j and r y < x<r, 

if i 7^ j and nj > x a i (4) 

if i = j. 



It is closely related to the contact (or adjacency) matrix, 
whose elements are unity for atom pairs within interac- 
tion range, and zero otherwise. Contact matrices have 
been used to describe the equilibrium structure of pro- 
teins in terms of amino acid contacts 15J, [l6| ; as described 



by Bahar et al. [17( and Sadoc [18( this idea may be ex- 



tended, through the dynamical matrix, to give a simple 
network model of vibrational motions. Here we suggest 
that the topology of the interactions in the chain around 
the critical point of collapse, described by T, may be 
linked to the dynamical effects that we have observed. 

The largest eigenvalue of T will be denoted as 7. Let 
a(A; 7) be the probability distribution of 7 at an interface 
A sampled by pathways started from A, and &(A; 7) the 
same quantity but sampled by pathways started from 
B. We found that a(A;7) is unimodal (approximately 
Gaussian) at all interfaces, with the mean value growing 
with A, and that the conformations at A with large 7 
are more likely to crystallise. The distribution b(X; 7) at 
interfaces far enough from the isocommittor is also uni- 
modal with similar properties, but becomes bimodal at 
the interfaces A close to the isocommittor. The critical 
value separating these two modes is denoted as j c . The 
insets of Fig. 2] show that pathways started from A, and 
reaching these energies, do not sample the population 
of low-eigenvalue states. The microscopic reversibility of 
our dynamics then implies that folding transition path- 
ways must cross X B _ 1 at 7 > j c . Indeed, the probability 
analysis in Fig. @] shows that pathways started at X B _ 1 
with 7 < 7 C have almost no chance to reach A. The dif- 
ference between the distributions a(A; 7) and 6(A;7) con- 
firms our suspicion that the forward-going and reverse- 
going ensembles of reactive trajectories are not identical 
(due to kinetic hindering) , and that the top eigenvalue of 
the dynamical or contact matrix may be able to discrim- 
inate between these families of trajectories. Incidentally, 
the equilibrium distribution of 7 at this energy, obtained 
by WL, is very similar to the nonreactive 7 < 7 C portion 
of Fig. HJb) . Similar observations apply to other eigen- 
values near the top of the spectrum of T. It is significant 
that typical 7 > 7 C configurations appear to have a com- 
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FIG. 4: Configurations at A^_i with largest eigenvalue 7 of 
the dynamical matrix lower than the critical value (7° = 11.9) 
have almost no chance to crystallize. Insets: (a) Unimodal 



distribution a(A ; 7) at the surface A 



-235 sampled by 



pathways started in A. (b) Bimodal distribution of b(\ B ;j) 
at the same surface (A^_i = —235) but sampled by pathways 
started in B. 



pact crystal nucleus with attached chains or loops, while 
for 7 < 7 C the same number of interactions are typically 
arranged in a single, less well ordered, cluster. 

It is worth mentioning that two structures on surfaces 
A close to the isocommittor have also been identified by 
Taylor et al. Q using the radius of gyration (R g ) as a sec- 
ond reaction coordinate. This result was also confirmed 
here. The correlation between R g and 7 was tested and 
found to be only weak. An analysis similar to that in 
Fig. [4] showed that 7 has significantly better predictive 
properties than R g . 

Why do we believe that the dynamical matrix, and 
the associated contact matrix, deserve further study? As 
mentioned above, they give a rather general connection 
between the topology defined by the interactions within 
a chain configuration and its dynamical evolution, in the 
approximation of an elastic network model. This has not 
only been used in the discussion_of pro teins to identify 
vibrational modes of oscillations 
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but also in the 
definition of nodes and an order parameter (a distance 
between nodes) i n dy namical network models of the fold- 
ing process itself [19( . The dynamical matrix is also being 
used for the description of glassy structures in colloidal 
systems [2(|. Most recently, the contact matrix of an 
atomic cluster has been used as a generator of order pa- 
rameters for metadynamics simulations 2l| . Our results 



clearly reinforce the view that the dynamical matrix is a 
simple object capturing successfully important topologi- 
cal or vibrational features of various interacting systems. 

To summarize, the transition of the homopolymer 
chain from the disordered globule to the crystal state 



has been simulated by dynamical forward flux sampling 
and brute force simulation. The results gave strong evi- 
dence that kinetic effects play an important role in the de- 
termination of the effective transition temperature. We 
then showed that the eigenvalues of the dynamical ma- 
trix yield further important information regarding the 
forward and reverse trajectories in the folding transition 
which complement the potential energy as an order pa- 
rameter. 
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